rm(list = ls())
setwd("~/Dropbox/Historical Displacement and Ethnicity - Afghanistan/REPOSITORY")

# Load necessary packages
library(ggplot2)
library(dplyr)
library(tidyr)
library(readr)

################################################################################
################################## Effective ###################################
################################################################################

# Read in the exported Stata results
results <- read.csv("replication/data/intermediate/FOGHORNresults_effective.csv")
results <- sapply(results, gsub, pattern = "^=", replacement = "")
colnames(results) <- results[1, ]
results <- results[-1, ]
results <- as.data.frame(results)
results[, 2:25] <- lapply(results[, 2:25], as.numeric)

# Mean of the residual category
means <- read.csv("replication/data/intermediate/FOGHORNresults_effective_means.csv")
means <- sapply(means, gsub, pattern = "^=", replacement = "")
colnames(means) <- means[1, ]
means <- means[-1, ]
means <- as.data.frame(t(means))
means[, 2:25] <- lapply(means[, 2:25], as.numeric)

# Extract coefficients and standard errors for each group
non_pashtun_displaced_effect     <- results[results$V1 == "nonpashtun_displaced", 2:25]
non_pashtun_displaced_se         <- results[which(results$V1 == "nonpashtun_displaced") + 1, 2:25]
non_pashtun_non_displaced_effect <- results[results$V1 == "nonpashtun_nondisplaced", 2:25]
non_pashtun_non_displaced_se     <- results[which(results$V1 == "nonpashtun_nondisplaced") + 1, 2:25]
displaced_pashtun_effect         <- results[results$V1 == "displaced_pashtun", 2:25]
displaced_pashtun_se             <- results[which(results$V1 == "displaced_pashtun") + 1, 2:25]
all_pashtun_north_effect         <- results[results$V1 == "all_pashtun_north", 2:25]
all_pashtun_north_se             <- results[which(results$V1 == "all_pashtun_north") + 1, 2:25]
displaced_pashtun2_effect        <- results[results$V1 == "displaced_pashtun2", 2:25]
displaced_pashtun2_se            <- results[which(results$V1 == "displaced_pashtun2") + 1, 2:25]

# Create a dataframe with the desired groups and their effects
effects_df <- data.frame(
  group  = c(
    rep("Displaced\nPashtuns", 24),
    rep("Non-Pashtuns\nin Displaced\nRegions", 24),
    rep("Non-Pashtuns\nNon-Displaced", 24),
    rep("All Pashtuns\nin the North", 24),
    rep("Displaced\nPashtuns\n(Alternative)", 24)
  ),
  effect = unlist(c(
    displaced_pashtun_effect,
    non_pashtun_displaced_effect,
    non_pashtun_non_displaced_effect,
    all_pashtun_north_effect,
    displaced_pashtun2_effect
  )),
  se     = unlist(c(
    displaced_pashtun_se,
    non_pashtun_displaced_se,
    non_pashtun_non_displaced_se,
    all_pashtun_north_se,
    displaced_pashtun2_se
  )),
  dv     = names(c(
    displaced_pashtun_effect,
    non_pashtun_displaced_effect,
    non_pashtun_non_displaced_effect,
    all_pashtun_north_effect,
    displaced_pashtun2_effect
  ))
)

effects_df1 <- effects_df[
  effects_df$dv %in% c(
    "prov_council1", 
    "cent_govt1", 
    "ulama1", 
    "locgovt_communityrelation1"
  ) & !is.na(effects_df$effect),
]

effects_df2 <- effects_df[
  effects_df$dv %in% c(
    "prov_council2", 
    "cent_govt2", 
    "ulama2", 
    "locgovt_communityrelation2"
  ) & !is.na(effects_df$effect),
]

##### plot effects_df1, without sigacts controls #####
effects_df1$dv <- factor(
  effects_df1$dv,
  levels = c(
    "prov_council1",
    "cent_govt1",
    "ulama1",
    "locgovt_communityrelation1"
  ),
  labels = c(
    paste0("Provincial\nCouncil\nEffective?\n(mean=", round(means$prov_council1, 3), ")"),
    paste0("Central\nGovernment\nEffective?\n(mean=", round(means$cent_govt1, 3), ")"),
    paste0("Ulama\nCouncils\nEffective?\n(mean=", round(means$ulama1, 3), ")"),
    paste0("Local\nGovernment\nRelations?\n(mean=", round(means$locgovt_communityrelation1, 3), ")")
  )
)

effects_df1 <- effects_df1[, c("dv", "group", "effect", "se")]
effects_df1 <- effects_df1[order(effects_df1$dv), ]

# Reordering the covariate and outcome categories
effects_df1$group <- factor(
  effects_df1$group,
  levels = c(
    "Non-Pashtuns\nNon-Displaced",
    "Non-Pashtuns\nin Displaced\nRegions",
    "Displaced\nPashtuns"
  )
)

pdf("replication/results/figures/FigureA15.pdf", width = 15, height = 5)
ggplot(effects_df1, aes(x = group, y = effect)) +
  geom_segment(
    aes(
      x      = group, 
      xend   = group, 
      y      = effect - 1.96 * se, 
      yend   = effect + 1.96 * se, 
      linetype = "95% Confidence Interval"
    ),
    linewidth = 2
  ) +
  #geom_segment(
  #  aes(
  #    x      = group, 
  #    xend   = group, 
  #    y      = effect - 1.645 * se, 
  #    yend   = effect + 1.645 * se, 
  #    linetype = "90% Confidence Interval"
  #  )
  #) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_wrap(~ dv, nrow = 1) +
  coord_flip() +
  xlab("") +
  ylab("") +  # "Compared to the baseline group, Non-Pashtuns Non-Displaced"
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_text(size = 16),
    strip.text         = element_text(size = 16),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom",
    panel.spacing.x    = unit(2, "lines")
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
dev.off()

##### plot effects_df2, with sigacts controls #####
effects_df2$dv <- factor(
  effects_df2$dv,
  levels = c(
    "prov_council2",
    "cent_govt2",
    "ulama2",
    "locgovt_communityrelation2"
  ),
  labels = c(
    paste0("Provincial\nCouncil\nEffective?\n(mean=", round(means$prov_council2, 3), ")"),
    paste0("Central\nGovernment\nEffective?\n(mean=", round(means$cent_govt2, 3), ")"),
    paste0("Ulama\nCouncils\nEffective?\n(mean=", round(means$ulama2, 3), ")"),
    paste0("Local\nGovernment\nRelations?\n(mean=", round(means$locgovt_communityrelation2, 3), ")")
  )
)

effects_df2 <- effects_df2[, c("dv", "group", "effect", "se")]
effects_df2 <- effects_df2[order(effects_df2$dv), ]

# Reordering the covariate and outcome categories
effects_df2$group <- factor(
  effects_df2$group,
  levels = c(
    "Non-Pashtuns\nNon-Displaced",
    "Non-Pashtuns\nin Displaced\nRegions",
    "Displaced\nPashtuns"
  )
)

pdf("replication/results/figures/Figure6.pdf", width = 15, height = 5)
ggplot(effects_df2, aes(x = group, y = effect)) +
  geom_segment(
    aes(
      x      = group,
      xend   = group,
      y      = effect - 1.96 * se,
      yend   = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    linewidth = 2
  ) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_wrap(~ dv, nrow = 1) +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_text(size = 16),
    strip.text         = element_text(size = 16),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom",
    panel.spacing.x    = unit(2, "lines")
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
dev.off()


###########################
### Alternative Measure ###
###########################
effects_df6 <- effects_df[
  effects_df$dv %in% c("prov_council6", "cent_govt6", "ulama6", "locgovt_communityrelation6") &
    !is.na(effects_df$effect),
]

effects_df6$dv <- factor(
  effects_df6$dv,
  levels = c("prov_council6", "cent_govt6", "ulama6", "locgovt_communityrelation6"),
  labels = c(
    paste0("Provincial\nCouncil\nEffective?\n(mean=", round(means$prov_council6, 3), ")"),
    paste0("Central\nGovernment\nEffective?\n(mean=", round(means$cent_govt6, 3), ")"),
    paste0("Ulama\nCouncils\nEffective?\n(mean=", round(means$ulama6, 3), ")"),
    paste0("Local\nGovernment\nRelations?\n(mean=", round(means$locgovt_communityrelation6, 3), ")")
  )
)

effects_df6 <- effects_df6[, c("dv", "group", "effect", "se")]
effects_df6 <- effects_df6[order(effects_df6$dv), ]

# Reordering the covariate and outcome categories
effects_df6$group <- factor(
  effects_df6$group,
  levels = c(
    "Non-Pashtuns\nNon-Displaced",
    "Non-Pashtuns\nin Displaced\nRegions",
    "Displaced\nPashtuns\n(Alternative)"
  )
)

pdf("replication/results/figures/FigureA11.pdf", width = 15, height = 5)
ggplot(effects_df6, aes(x = group, y = effect)) +
  geom_segment(
    aes(
      x        = group,
      xend     = group,
      y        = effect - 1.96 * se,
      yend     = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    linewidth = 2
  ) +
  # geom_segment(aes(x = group, xend = group,
  #                  y = effect - 1.645 * se,
  #                  yend = effect + 1.645 * se,
  #                  linetype = "90% Confidence Interval")) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_wrap(~ dv, nrow = 1) +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_text(size = 16),
    strip.text         = element_text(size = 16),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom",
    panel.spacing.x    = unit(2, "lines")
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
dev.off()

###################################
### different comparison groups ###
###################################
means <- as.data.frame(
  cbind(colnames(means)[2:25], t(means[, 2:25]))
)
colnames(means) <- c("dv", "mean")
means$mean <- as.numeric(means$mean)

effects_df <- effects_df[
  !is.na(effects_df$effect) &
    substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) %in% c(2, 3, 4, 5) &
    !effects_df$group %in% c("Non-Pashtuns\nNon-Displaced", "Non-Pashtuns\nin Displaced\nRegions"),
]

effects_df <- merge(effects_df, means, by = "dv", all.x = TRUE)

# specify the comparison group (explanatory variables)
effects_df$group[
  substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) == 2
] <- "As presented in\nthe main manuscript"

effects_df$group[
  substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) == 3
] <- "Displaced Pashtuns\nvs.\nPashtuns in the South"

effects_df$group[
  substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) == 4
] <- "Pashtuns in the North\nvs.\nAll Others"

effects_df$group[
  substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) == 5
] <- "Displaced Pashtuns\nvs.\nOthers in the North,\noutside the Displaced Region"

effects_df$group <- factor(
  effects_df$group,
  levels = c(
    "As presented in\nthe main manuscript",
    "Displaced Pashtuns\nvs.\nPashtuns in the South",
    "Pashtuns in the North\nvs.\nAll Others",
    "Displaced Pashtuns\nvs.\nOthers in the North,\noutside the Displaced Region"
  )
)

effects_df <- effects_df[order(effects_df$group, decreasing = TRUE), ]

# (dependent variables)
effects_df$dv <- ifelse(
  effects_df$dv %in% c("prov_council2", "prov_council3", "prov_council4", "prov_council5"),
  "Provincial\nCouncil\nEffective?",
  effects_df$dv
)
effects_df$dv <- ifelse(
  effects_df$dv %in% c("cent_govt2", "cent_govt3", "cent_govt4", "cent_govt5"),
  "Central\nGovernment\nEffective?",
  effects_df$dv
)
effects_df$dv <- ifelse(
  effects_df$dv %in% c("ulama2", "ulama3", "ulama4", "ulama5"),
  "Ulama\nCouncils\nEffective?",
  effects_df$dv
)
effects_df$dv <- ifelse(
  effects_df$dv %in% c("locgovt_communityrelation2", "locgovt_communityrelation3", "locgovt_communityrelation4", "locgovt_communityrelation5"),
  "Local\nGovernment\nRelations?",
  effects_df$dv
)

effects_df$dv <- factor(
  effects_df$dv,
  levels = c(
    "Provincial\nCouncil\nEffective?",
    "Central\nGovernment\nEffective?",
    "Ulama\nCouncils\nEffective?",
    "Local\nGovernment\nRelations?"
  )
)

pdf("replication/results/figures/FigureA7.pdf", width = 16, height = 12)
ggplot(effects_df, aes(x = 1, y = effect)) +
  geom_segment(
    aes(
      xend     = 1,
      y        = effect - 1.96 * se,
      yend     = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    linewidth = 2
  ) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_grid(group ~ dv, switch = "y") +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_blank(),  # Remove y-axis text
    axis.ticks.y       = element_blank(),  # Remove y-axis ticks
    strip.text         = element_text(size = 16),
    strip.text.x       = element_text(size = 16, margin = margin(t = 20, b = 20)),
    strip.text.y.left  = element_text(size = 16, angle = 0, hjust = 0.5),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom",
    panel.spacing      = unit(2, "lines")
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1)) +
  geom_text(
    data    = effects_df,
    aes(
      x     = 1,
      y     = Inf,
      label = paste0("mean=", round(mean, 3))
    ),
    hjust       = 1.05,
    vjust       = 4,
    size        = 5,
    inherit.aes = FALSE
  )
dev.off()

################################################################################
################################## Legitimacy ##################################
################################################################################
rm(list = ls())

# Read in the exported Stata results
results <- read.csv("replication/data/intermediate/FOGHORNresults_legit.csv")
results <- sapply(results, gsub, pattern = "^=", replacement = "")
colnames(results) <- results[1, ]
results <- results[-1, ]
results <- as.data.frame(results)
results[, 2:25] <- lapply(results[, 2:25], as.numeric)

# Mean of the residual category
means <- read.csv("replication/data/intermediate/FOGHORNresults_legit_means.csv")
means <- sapply(means, gsub, pattern = "^=", replacement = "")
colnames(means) <- means[1, ]
means <- means[-1, ]
means <- as.data.frame(t(means))
means[, 2:25] <- lapply(means[, 2:25], as.numeric)

# Extract coefficients and standard errors for each group
non_pashtun_displaced_effect     <- results[results$V1 == "nonpashtun_displaced", 2:25]
non_pashtun_displaced_se         <- results[which(results$V1 == "nonpashtun_displaced") + 1, 2:25]
non_pashtun_non_displaced_effect <- results[results$V1 == "nonpashtun_nondisplaced", 2:25]
non_pashtun_non_displaced_se     <- results[which(results$V1 == "nonpashtun_nondisplaced") + 1, 2:25]
displaced_pashtun_effect         <- results[results$V1 == "displaced_pashtun", 2:25]
displaced_pashtun_se             <- results[which(results$V1 == "displaced_pashtun") + 1, 2:25]
all_pashtun_north_effect         <- results[results$V1 == "all_pashtun_north", 2:25]
all_pashtun_north_se             <- results[which(results$V1 == "all_pashtun_north") + 1, 2:25]
displaced_pashtun2_effect        <- results[results$V1 == "displaced_pashtun2", 2:25]
displaced_pashtun2_se            <- results[which(results$V1 == "displaced_pashtun2") + 1, 2:25]

# Create a dataframe with the desired groups and their effects
effects_df <- data.frame(
  group  = c(
    rep("Displaced\nPashtuns", 24),
    rep("Non-Pashtuns\nin Displaced\nRegions", 24),
    rep("Non-Pashtuns\nNon-Displaced", 24),
    rep("All Pashtuns\nin the North", 24),
    rep("Displaced\nPashtuns\n(Alternative)", 24)
  ),
  effect = unlist(c(
    displaced_pashtun_effect,
    non_pashtun_displaced_effect,
    non_pashtun_non_displaced_effect,
    all_pashtun_north_effect,
    displaced_pashtun2_effect
  )),
  se     = unlist(c(
    displaced_pashtun_se,
    non_pashtun_displaced_se,
    non_pashtun_non_displaced_se,
    all_pashtun_north_se,
    displaced_pashtun2_se
  )),
  dv     = names(c(
    displaced_pashtun_effect,
    non_pashtun_displaced_effect,
    non_pashtun_non_displaced_effect,
    all_pashtun_north_effect,
    displaced_pashtun2_effect
  ))
)

effects_df1 <- effects_df[
  effects_df$dv %in% c("provgovt1", "jirga1", "centgovt1", "judicial1") &
    !is.na(effects_df$effect),
]
effects_df2 <- effects_df[
  effects_df$dv %in% c("provgovt2", "jirga2", "centgovt2", "judicial2") &
    !is.na(effects_df$effect),
]

##### plot effects_df1, without sigacts controls #####
effects_df1$dv <- factor(
  effects_df1$dv,
  levels = c("provgovt1", "jirga1", "centgovt1", "judicial1"),
  labels = c(
    paste0("Provincial\nGovernment\nLegitimate?\n(mean=", round(means$provgovt1, 3), ")"),
    paste0("Wolesi\nJirga\nLegitimate?\n(mean=", round(means$jirga1, 3), ")"),
    paste0("Central\nGovernment\nLegitimate?\n(mean=", round(means$centgovt1, 3), ")"),
    paste0("Government\nCourts\nFair?\n(mean=", round(means$judicial1, 3), ")")
  )
)

effects_df1 <- effects_df1[, c("dv", "group", "effect", "se")]
effects_df1 <- effects_df1[order(effects_df1$dv), ]

# Reordering the covariate and outcome categories
effects_df1$group <- factor(
  effects_df1$group,
  levels = c(
    "Non-Pashtuns\nNon-Displaced",
    "Non-Pashtuns\nin Displaced\nRegions",
    "Displaced\nPashtuns"
  )
)

pdf("replication/results/figures/FigureA16.pdf", width = 15, height = 5)
ggplot(effects_df1, aes(x = group, y = effect)) +
  geom_segment(
    aes(
      x        = group,
      xend     = group,
      y        = effect - 1.96 * se,
      yend     = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    linewidth = 2
  ) +
  # geom_segment(aes(x = group, xend = group,
  #                  y = effect - 1.645 * se,
  #                  yend = effect + 1.645 * se,
  #                  linetype = "90% Confidence Interval")) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_wrap(~dv, nrow = 1) +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_text(size = 16),
    strip.text         = element_text(size = 16),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom",
    panel.spacing.x    = unit(2, "lines")
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
dev.off()

##### plot effects_df2, with sigacts controls #####
effects_df2$dv <- factor(
  effects_df2$dv,
  levels = c("provgovt2", "jirga2", "centgovt2", "judicial2"),
  labels = c(
    paste0("Provincial\nGovernment\nLegitimate?\n(mean=", round(means$provgovt2, 3), ")"),
    paste0("Wolesi\nJirga\nLegitimate?\n(mean=", round(means$jirga2, 3), ")"),
    paste0("Central\nGovernment\nLegitimate?\n(mean=", round(means$centgovt2, 3), ")"),
    paste0("Government\nCourts\nFair?\n(mean=", round(means$judicial2, 3), ")")
  )
)

effects_df2 <- effects_df2[, c("dv", "group", "effect", "se")]
effects_df2 <- effects_df2[order(effects_df2$dv), ]

# Reordering the covariate and outcome categories
effects_df2$group <- factor(
  effects_df2$group,
  levels = c(
    "Non-Pashtuns\nNon-Displaced",
    "Non-Pashtuns\nin Displaced\nRegions",
    "Displaced\nPashtuns"
  )
)

pdf("replication/results/figures/Figure7.pdf", width = 15, height = 5)
ggplot(effects_df2, aes(x = group, y = effect)) +
  geom_segment(
    aes(
      x        = group,
      xend     = group,
      y        = effect - 1.96 * se,
      yend     = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    linewidth = 2
  ) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_wrap(~dv, nrow = 1) +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_text(size = 16),
    strip.text         = element_text(size = 16),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom",
    panel.spacing.x    = unit(2, "lines")
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
dev.off()


###########################
### Alternative Measure ###
###########################
effects_df6 <- effects_df[
  effects_df$dv %in% c("provgovt6", "jirga6", "centgovt6", "judicial6") &
    !is.na(effects_df$effect),
]

effects_df6$dv <- factor(
  effects_df6$dv,
  levels = c("provgovt6", "jirga6", "centgovt6", "judicial6"),
  labels = c(
    paste0("Provincial\nGovernment\nLegitimate?\n(mean=", round(means$provgovt6, 3), ")"),
    paste0("Wolesi\nJirga\nLegitimate?\n(mean=", round(means$jirga6, 3), ")"),
    paste0("Central\nGovernment\nLegitimate?\n(mean=", round(means$centgovt6, 3), ")"),
    paste0("Government\nCourts\nFair?\n(mean=", round(means$judicial6, 3), ")")
  )
)

effects_df6 <- effects_df6[, c("dv", "group", "effect", "se")]
effects_df6 <- effects_df6[order(effects_df6$dv), ]

# Reordering the covariate and outcome categories
effects_df6$group <- factor(
  effects_df6$group,
  levels = c(
    "Non-Pashtuns\nNon-Displaced",
    "Non-Pashtuns\nin Displaced\nRegions",
    "Displaced\nPashtuns\n(Alternative)"
  )
)

pdf("replication/results/figures/FigureA12.pdf", width = 15, height = 5)
ggplot(effects_df6, aes(x = group, y = effect)) +
  geom_segment(
    aes(
      x        = group,
      xend     = group,
      y        = effect - 1.96 * se,
      yend     = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    linewidth = 2
  ) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_wrap(~ dv, nrow = 1) +
  coord_flip() +
  xlab("") +
  ylab("") +  # "Compared to the baseline group, Non-Pashtuns Non-Displaced"
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_text(size = 16),
    strip.text         = element_text(size = 16),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom",
    panel.spacing.x    = unit(2, "lines")
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
dev.off()


###################################
### different comparison groups ###
###################################
means <- as.data.frame(
  cbind(colnames(means)[2:25], t(means[, 2:25]))
)
colnames(means) <- c("dv", "mean")
means$mean <- as.numeric(means$mean)

effects_df <- effects_df[
  !is.na(effects_df$effect) &
    substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) %in% c(2, 3, 4, 5) &
    !effects_df$group %in% c(
      "Non-Pashtuns\nNon-Displaced",
      "Non-Pashtuns\nin Displaced\nRegions"
    ),
]

effects_df <- merge(effects_df, means, by = "dv", all.x = TRUE)

# specify the comparison group (explanatory variables)
effects_df$group[
  substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) == 2
] <- "As presented in\nthe main manuscript"
effects_df$group[
  substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) == 3
] <- "Displaced Pashtuns\nvs.\nPashtuns in the South"
effects_df$group[
  substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) == 4
] <- "Pashtuns in the North\nvs.\nAll Others"
effects_df$group[
  substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) == 5
] <- "Displaced Pashtuns\nvs.\nOthers in the North,\noutside the Displaced Region"

effects_df$group <- factor(
  effects_df$group,
  levels = c(
    "As presented in\nthe main manuscript",
    "Displaced Pashtuns\nvs.\nPashtuns in the South",
    "Pashtuns in the North\nvs.\nAll Others",
    "Displaced Pashtuns\nvs.\nOthers in the North,\noutside the Displaced Region"
  )
)
effects_df <- effects_df[order(effects_df$group, decreasing = TRUE), ]

# (dependent variables)
effects_df$dv <- ifelse(
  effects_df$dv %in% c("provgovt2", "provgovt3", "provgovt4", "provgovt5"),
  "Provincial\nGovernment\nLegitimate?",
  effects_df$dv
)
effects_df$dv <- ifelse(
  effects_df$dv %in% c("jirga2", "jirga3", "jirga4", "jirga5"),
  "Wolesi\nJirga\nLegitimate?",
  effects_df$dv
)
effects_df$dv <- ifelse(
  effects_df$dv %in% c("centgovt2", "centgovt3", "centgovt4", "centgovt5"),
  "Central\nGovernment\nLegitimate?",
  effects_df$dv
)
effects_df$dv <- ifelse(
  effects_df$dv %in% c("judicial2", "judicial3", "judicial4", "judicial5"),
  "Government\nCourts\nFair?",
  effects_df$dv
)

effects_df$dv <- factor(
  effects_df$dv,
  levels = c(
    "Provincial\nGovernment\nLegitimate?",
    "Wolesi\nJirga\nLegitimate?",
    "Central\nGovernment\nLegitimate?",
    "Government\nCourts\nFair?"
  )
)

pdf("replication/results/figures/FigureA8.pdf", width = 16, height = 12)
ggplot(effects_df, aes(x = 1, y = effect)) +
  geom_segment(
    aes(
      xend     = 1,
      y        = effect - 1.96 * se,
      yend     = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    linewidth = 2
  ) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_grid(group ~ dv, switch = "y") +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_blank(),  # Remove y-axis text
    axis.ticks.y       = element_blank(),  # Remove y-axis ticks
    strip.text         = element_text(size = 16),
    strip.text.x       = element_text(size = 16, margin = margin(t = 20, b = 20)),
    strip.text.y.left  = element_text(size = 16, angle = 0, hjust = 0.5),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom",
    panel.spacing      = unit(2, "lines")
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1)) +
  geom_text(
    data    = effects_df,
    aes(
      x     = 1,
      y     = Inf,
      label = paste0("mean=", round(mean, 3))
    ),
    hjust       = 1.05,
    vjust       = 4,
    size        = 5,
    inherit.aes = FALSE
  )
dev.off()

################################################################################
################################### Taliban ####################################
################################################################################
rm(list = ls())

# Read in the exported Stata results
results <- read.csv("replication/data/intermediate/FOGHORNresults_taliban.csv")
results <- sapply(results, gsub, pattern = "^=", replacement = "")
colnames(results) <- results[1, ]
results <- results[-1, ]
results <- as.data.frame(results)
results[, 2:25] <- lapply(results[, 2:25], as.numeric)

# Mean of the residual category
means <- read.csv("replication/data/intermediate/FOGHORNresults_taliban_means.csv")
means <- sapply(means, gsub, pattern = "^=", replacement = "")
colnames(means) <- means[1, ]
means <- means[-1, ]
means <- as.data.frame(t(means))
means[, 2:25] <- lapply(means[, 2:25], as.numeric)

# Extract coefficients and standard errors for each group
non_pashtun_displaced_effect     <- results[results$V1 == "nonpashtun_displaced", 2:25]
non_pashtun_displaced_se         <- results[which(results$V1 == "nonpashtun_displaced") + 1, 2:25]
non_pashtun_non_displaced_effect <- results[results$V1 == "nonpashtun_nondisplaced", 2:25]
non_pashtun_non_displaced_se     <- results[which(results$V1 == "nonpashtun_nondisplaced") + 1, 2:25]
displaced_pashtun_effect         <- results[results$V1 == "displaced_pashtun", 2:25]
displaced_pashtun_se             <- results[which(results$V1 == "displaced_pashtun") + 1, 2:25]
all_pashtun_north_effect         <- results[results$V1 == "all_pashtun_north", 2:25]
all_pashtun_north_se             <- results[which(results$V1 == "all_pashtun_north") + 1, 2:25]
displaced_pashtun2_effect        <- results[results$V1 == "displaced_pashtun2", 2:25]
displaced_pashtun2_se            <- results[which(results$V1 == "displaced_pashtun2") + 1, 2:25]

# Create a dataframe with the desired groups and their effects
effects_df <- data.frame(
  group  = c(
    rep("Displaced\nPashtuns", 24),
    rep("Non-Pashtuns\nin Displaced\nRegions", 24),
    rep("Non-Pashtuns\nNon-Displaced", 24),
    rep("All Pashtuns\nin the North", 24),
    rep("Displaced\nPashtuns\n(Alternative)", 24)
  ),
  effect = unlist(c(
    displaced_pashtun_effect,
    non_pashtun_displaced_effect,
    non_pashtun_non_displaced_effect,
    all_pashtun_north_effect,
    displaced_pashtun2_effect
  )),
  se     = unlist(c(
    displaced_pashtun_se,
    non_pashtun_displaced_se,
    non_pashtun_non_displaced_se,
    all_pashtun_north_se,
    displaced_pashtun2_se
  )),
  dv     = names(c(
    displaced_pashtun_effect,
    non_pashtun_displaced_effect,
    non_pashtun_non_displaced_effect,
    all_pashtun_north_effect,
    displaced_pashtun2_effect
  ))
)

effects_df1 <- effects_df[
  effects_df$dv %in% c("more_influence1", "communityrelation1", "trust1", "courts_fair1") &
    !is.na(effects_df$effect),
]
effects_df2 <- effects_df[
  effects_df$dv %in% c("more_influence2", "communityrelation2", "trust2", "courts_fair2") &
    !is.na(effects_df$effect),
]

##### plot effects_df1, without sigacts controls #####
effects_df1$dv <- factor(
  effects_df1$dv,
  levels = c("more_influence1", "communityrelation1", "trust1", "courts_fair1"),
  labels = c(
    paste0("Does the Taliban\nHave More\nInfluence than\nthe Government?\n(mean=", round(means$more_influence1, 3), ")"),
    paste0("Are Taliban\nRelations with\nYour Community\nFavorable?\n(mean=", round(means$communityrelation1, 3), ")"),
    paste0("Do You Trust\nStatements\nMade by\nthe Taliban?\n(mean=", round(means$trust1, 3), ")"),
    paste0("Are\nTaliban Courts\nFair?\n(mean=", round(means$courts_fair1, 3), ")")
  )
)

effects_df1 <- effects_df1[, c("dv", "group", "effect", "se")]
effects_df1 <- effects_df1[order(effects_df1$dv), ]

# Reordering the covariate and outcome categories
effects_df1$group <- factor(
  effects_df1$group,
  levels = c(
    "Non-Pashtuns\nNon-Displaced",
    "Non-Pashtuns\nin Displaced\nRegions",
    "Displaced\nPashtuns"
  )
)

pdf("replication/results/figures/FigureA17.pdf", width = 15, height = 5)
ggplot(effects_df1, aes(x = group, y = effect)) +
  geom_segment(
    aes(
      x        = group,
      xend     = group,
      y        = effect - 1.96 * se,
      yend     = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    linewidth = 2
  ) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_wrap(~dv, nrow = 1) +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_text(size = 16),
    strip.text         = element_text(size = 16),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom",
    panel.spacing.x    = unit(2, "lines")
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
dev.off()

##### plot effects_df2, with sigacts controls #####
effects_df2$dv <- factor(
  effects_df2$dv,
  levels = c("more_influence2", "communityrelation2", "trust2", "courts_fair2"),
  labels = c(
    paste0("Does the Taliban\nHave More\nInfluence than\nthe Government?\n(mean=", round(means$more_influence2, 3), ")"),
    paste0("Are Taliban\nRelations with\nYour Community\nFavorable?\n(mean=", round(means$communityrelation2, 3), ")"),
    paste0("Do You Trust\nStatements\nMade by\nthe Taliban?\n(mean=", round(means$trust2, 3), ")"),
    paste0("Are\nTaliban Courts\nFair?\n(mean=", round(means$courts_fair2, 3), ")")
  )
)

effects_df2 <- effects_df2[, c("dv", "group", "effect", "se")]
effects_df2 <- effects_df2[order(effects_df2$dv), ]

# Reordering the covariate and outcome categories
effects_df2$group <- factor(
  effects_df2$group,
  levels = c(
    "Non-Pashtuns\nNon-Displaced",
    "Non-Pashtuns\nin Displaced\nRegions",
    "Displaced\nPashtuns"
  )
)

pdf("replication/results/figures/Figure8.pdf", width = 15, height = 5)
ggplot(effects_df2, aes(x = group, y = effect)) +
  geom_segment(
    aes(
      x        = group,
      xend     = group,
      y        = effect - 1.96 * se,
      yend     = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    linewidth = 2
  ) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_wrap(~dv, nrow = 1) +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_text(size = 16),
    strip.text         = element_text(size = 16),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom",
    panel.spacing.x    = unit(2, "lines")
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
dev.off()

###########################
### Alternative Measure ###
###########################
effects_df6 <- effects_df[
  effects_df$dv %in% c("more_influence6", "communityrelation6", "trust6", "courts_fair6") &
    !is.na(effects_df$effect),
]

effects_df6$dv <- factor(
  effects_df6$dv,
  levels = c("more_influence6", "communityrelation6", "trust6", "courts_fair6"),
  labels = c(
    paste0("Does the Taliban\nHave More\nInfluence than\nthe Government?\n(mean=", round(means$more_influence6, 3), ")"),
    paste0("Are Taliban\nRelations with\nYour Community\nFavorable?\n(mean=", round(means$communityrelation6, 3), ")"),
    paste0("Do You Trust\nStatements\nMade by\nthe Taliban?\n(mean=", round(means$trust6, 3), ")"),
    paste0("Are\nTaliban Courts\nFair?\n(mean=", round(means$courts_fair6, 3), ")")
  )
)

effects_df6 <- effects_df6[, c("dv", "group", "effect", "se")]
effects_df6 <- effects_df6[order(effects_df6$dv), ]

# Reordering the covariate and outcome categories
effects_df6$group <- factor(
  effects_df6$group,
  levels = c(
    "Non-Pashtuns\nNon-Displaced",
    "Non-Pashtuns\nin Displaced\nRegions",
    "Displaced\nPashtuns\n(Alternative)"
  )
)

pdf("replication/results/figures/FigureA13.pdf", width = 15, height = 5)
ggplot(effects_df6, aes(x = group, y = effect)) +
  geom_segment(
    aes(
      x        = group,
      xend     = group,
      y        = effect - 1.96 * se,
      yend     = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    linewidth = 2
  ) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_wrap(~dv, nrow = 1) +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_text(size = 16),
    strip.text         = element_text(size = 16),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom",
    panel.spacing.x    = unit(2, "lines")
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
dev.off()


###################################
### different comparison groups ###
###################################
means <- as.data.frame(
  cbind(colnames(means)[2:25], t(means[, 2:25]))
)
colnames(means) <- c("dv", "mean")
means$mean <- as.numeric(means$mean)

effects_df <- effects_df[
  !is.na(effects_df$effect) &
    substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) %in% c(2, 3, 4, 5) &
    !effects_df$group %in% c(
      "Non-Pashtuns\nNon-Displaced",
      "Non-Pashtuns\nin Displaced\nRegions"
    ),
]

effects_df <- merge(effects_df, means, by = "dv", all.x = TRUE)

# specify the comparison group (explanatory variables)
effects_df$group[
  substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) == 2
] <- "As presented in\nthe main manuscript"

effects_df$group[
  substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) == 3
] <- "Displaced Pashtuns\nvs.\nPashtuns in the South"

effects_df$group[
  substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) == 4
] <- "Pashtuns in the North\nvs.\nAll Others"

effects_df$group[
  substr(effects_df$dv, nchar(effects_df$dv), nchar(effects_df$dv)) == 5
] <- "Displaced Pashtuns\nvs.\nOthers in the North,\noutside the Displaced Region"

effects_df$group <- factor(
  effects_df$group,
  levels = c(
    "As presented in\nthe main manuscript",
    "Displaced Pashtuns\nvs.\nPashtuns in the South",
    "Pashtuns in the North\nvs.\nAll Others",
    "Displaced Pashtuns\nvs.\nOthers in the North,\noutside the Displaced Region"
  )
)

effects_df <- effects_df[order(effects_df$group, decreasing = TRUE), ]

# (dependent variables)
effects_df$dv <- ifelse(
  effects_df$dv %in% c("more_influence2", "more_influence3", "more_influence4", "more_influence5"),
  "Does the Taliban\nHave More\nInfluence than\nthe Government?",
  effects_df$dv
)
effects_df$dv <- ifelse(
  effects_df$dv %in% c("communityrelation2", "communityrelation3", "communityrelation4", "communityrelation5"),
  "Are Taliban\nRelations with\nYour Community\nFavorable?",
  effects_df$dv
)
effects_df$dv <- ifelse(
  effects_df$dv %in% c("trust2", "trust3", "trust4", "trust5"),
  "Do You Trust\nStatements\nMade by\nthe Taliban?",
  effects_df$dv
)
effects_df$dv <- ifelse(
  effects_df$dv %in% c("courts_fair2", "courts_fair3", "courts_fair4", "courts_fair5"),
  "Are\nTaliban Courts\nFair?",
  effects_df$dv
)

effects_df$dv <- factor(
  effects_df$dv,
  levels = c(
    "Does the Taliban\nHave More\nInfluence than\nthe Government?",
    "Are Taliban\nRelations with\nYour Community\nFavorable?",
    "Do You Trust\nStatements\nMade by\nthe Taliban?",
    "Are\nTaliban Courts\nFair?"
  )
)

pdf("replication/results/figures/FigureA9.pdf", width = 16, height = 12)
ggplot(effects_df, aes(x = 1, y = effect)) +
  geom_segment(
    aes(
      xend     = 1,
      y        = effect - 1.96 * se,
      yend     = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    linewidth = 2
  ) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_grid(group ~ dv, switch = "y") +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_blank(),  # Remove y-axis text
    axis.ticks.y       = element_blank(),  # Remove y-axis ticks
    strip.text         = element_text(size = 16),
    strip.text.x       = element_text(size = 16, margin = margin(t = 20, b = 20)),
    strip.text.y.left  = element_text(size = 16, angle = 0, hjust = 0.5),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom",
    panel.spacing      = unit(2, "lines")
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1)) +
  geom_text(
    data    = effects_df,
    aes(
      x     = 1,
      y     = Inf,
      label = paste0("mean=", round(mean, 3))
    ),
    hjust   = 1.05,
    vjust   = 4,
    size    = 5,
    inherit.aes = FALSE
  )
dev.off()

################################################################################
##################################### Fled Home #################################
################################################################################
rm(list = ls())

# Read in the exported Stata results
results <- read.csv("replication/data/intermediate/FOGHORNresults_fledhome.csv")
results <- sapply(results, gsub, pattern = "^=", replacement = "")
colnames(results) <- results[1, ]
results <- results[-1, ]
results <- as.data.frame(results)
results[, 2:13] <- lapply(results[, 2:13], as.numeric)

# Mean of the residual category
means <- read.csv("replication/data/intermediate/FOGHORNresults_fledhome_means.csv")
means <- sapply(means, gsub, pattern = "^=", replacement = "")
colnames(means) <- means[1, ]
means <- means[-1, ]
means <- as.data.frame(t(means))
means[, 2:13] <- lapply(means[, 2:13], as.numeric)

# Extract coefficients and standard errors for each group
non_pashtun_displaced_effect     <- results[results$V1 == "nonpashtun_displaced", 2:13]
non_pashtun_displaced_se         <- results[which(results$V1 == "nonpashtun_displaced") + 1, 2:13]
non_pashtun_non_displaced_effect <- results[results$V1 == "nonpashtun_nondisplaced", 2:13]
non_pashtun_non_displaced_se     <- results[which(results$V1 == "nonpashtun_nondisplaced") + 1, 2:13]
displaced_pashtun_effect         <- results[results$V1 == "displaced_pashtun", 2:13]
displaced_pashtun_se             <- results[which(results$V1 == "displaced_pashtun") + 1, 2:13]

# Create a dataframe with the desired groups and their effects
effects_df <- data.frame(
  group  = c(
    rep("Displaced\nPashtuns",            12),
    rep("Non-Pashtuns\nin Displaced\nRegions", 12),
    rep("Non-Pashtuns\nNon-Displaced",    12)
  ),
  effect = unlist(c(
    displaced_pashtun_effect,
    non_pashtun_displaced_effect,
    non_pashtun_non_displaced_effect
  )),
  se     = unlist(c(
    displaced_pashtun_se,
    non_pashtun_displaced_se,
    non_pashtun_non_displaced_se
  )),
  dv     = names(c(
    displaced_pashtun_effect,
    non_pashtun_displaced_effect,
    non_pashtun_non_displaced_effect
  ))
)

library(stringr)

effects_df2 <- effects_df %>%
  filter(grepl("polviol2|jobs2|anyreason2|community_safety2", dv)) %>%
  mutate(
    agegroup = case_when(
      str_ends(dv, "allage") ~ "All Ages",
      str_ends(dv, "old")    ~ "Older than 50",
      str_ends(dv, "young")  ~ "Younger than 50"
    ),
    dv_main = str_replace(dv, "(allage|old|young)$", "")
  )

# Factor dv_main with desired labels
effects_df2$dv_main <- factor(
  effects_df2$dv_main,
  levels = c("polviol2", "jobs2", "anyreason2", "community_safety2"),
  labels = c(
    paste0(
      "Fled Home\nDue to\nPolitical Violence\n",
      "(mean, all ages=", round(means$polviol2allage, 3), ")\n",
      "(mean, '>=50'=",        round(means$polviol2old,   3), ")\n",
      "(mean, '<50'=",         round(means$polviol2young, 3), ")"
    ),
    paste0(
      "Fled Home\nDue to\nEconomic Reasons\n",
      "(mean, all ages=", round(means$jobs2allage, 3), ")\n",
      "(mean, '>=50'=",        round(means$jobs2old,   3), ")\n",
      "(mean, '<50'=",         round(means$jobs2young, 3), ")"
    ),
    paste0(
      "Moved for\nAny Reason\n",
      "(mean, all ages=", round(means$anyreason2allage, 3), ")\n",
      "(mean, '>=50'=",        round(means$anyreason2old,   3), ")\n",
      "(mean, '<50'=",         round(means$anyreason2young, 3), ")"
    ),
    paste0(
      "Feel Safe in\nCommunity\n",
      "(mean, all ages=", round(means$community_safety2allage, 3), ")\n",
      "(mean, '>=50'=",        round(means$community_safety2old,   3), ")\n",
      "(mean, '<50'=",         round(means$community_safety2young, 3), ")"
    )
  )
)

effects_df2 <- effects_df2[, c("dv_main", "agegroup", "group", "effect", "se")]
effects_df2$group <- factor(
  effects_df2$group,
  levels = c(
    "Non-Pashtuns\nNon-Displaced",
    "Non-Pashtuns\nin Displaced\nRegions",
    "Displaced\nPashtuns"
  )
)

pdf("replication/results/figures/Figure5.pdf", width = 12, height = 6)
ggplot(effects_df2, aes(x = group, y = effect, color = agegroup, linetype = agegroup)) +
  geom_hline(yintercept = 0, color = "red4") +
  
  # All Ages (top line, no nudge)
  geom_segment(
    data     = subset(effects_df2, agegroup == "All Ages"),
    aes(xend = group, y = effect - 1.96 * se, yend = effect + 1.96 * se),
    linewidth = 1,
    position  = position_nudge(x = 0)
  ) +
  
  # Younger than 50 (below All Ages)
  geom_segment(
    data     = subset(effects_df2, agegroup == "Younger than 50"),
    aes(xend = group, y = effect - 1.96 * se, yend = effect + 1.96 * se),
    linewidth = 1,
    position  = position_nudge(x = -0.2)
  ) +
  
  # Older than 50 (below Younger than 50)
  geom_segment(
    data     = subset(effects_df2, agegroup == "Older than 50"),
    aes(xend = group, y = effect - 1.96 * se, yend = effect + 1.96 * se),
    linewidth = 1,
    position  = position_nudge(x = -0.4)
  ) +
  
  facet_wrap(~ dv_main, nrow = 1) +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_text(size = 16),
    strip.text         = element_text(size = 16),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_text(size = 16),
    legend.text        = element_text(size = 16),
    legend.position    = "bottom"
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1)) +
  
  # Manually set colors and linetypes for easy grayscale differentiation
  scale_color_manual(
    name   = "95% Confidence Interval",
    values = c(
      "All Ages"        = "black",
      "Younger than 50" = "green",
      "Older than 50"   = "blue"
    )
  ) +
  scale_linetype_manual(
    name   = "95% Confidence Interval",
    values = c(
      "All Ages"        = "solid",
      "Younger than 50" = "dotted",
      "Older than 50"   = "dashed"
    )
  )
dev.off()


################################################################################
############################### Economic Status ################################
################################################################################
rm(list = ls())

# Read in the exported Stata results
results <- read.csv("replication/data/intermediate/FOGHORNresults_econstatus.csv")
results <- sapply(results, gsub, pattern = "^=", replacement = "")
colnames(results) <- results[1, ]
results <- results[-1, ]
results <- as.data.frame(results)
results[, 2:3] <- lapply(results[, 2:3], as.numeric)

# Extract coefficients and standard errors for each group
non_pashtun_displaced_effect     <- results[results$V1 == "nonpashtun_displaced", 3]
non_pashtun_displaced_se         <- results[which(results$V1 == "nonpashtun_displaced") + 1, 3]
non_pashtun_non_displaced_effect <- results[results$V1 == "nonpashtun_nondisplaced", 3]
non_pashtun_non_displaced_se     <- results[which(results$V1 == "nonpashtun_nondisplaced") + 1, 3]
displaced_pashtun_effect         <- results[results$V1 == "displaced_pashtun", 3]
displaced_pashtun_se             <- results[which(results$V1 == "displaced_pashtun") + 1, 3]

# Create a dataframe with the desired groups and their effects
effects_df <- data.frame(
  group  = c(
    rep("Displaced\nPashtuns",            1),
    rep("Non-Pashtuns\nin Displaced\nRegions", 1),
    rep("Non-Pashtuns\nNon-Displaced",    1)
  ),
  effect = unlist(c(
    displaced_pashtun_effect,
    non_pashtun_displaced_effect,
    non_pashtun_non_displaced_effect
  )),
  se     = unlist(c(
    displaced_pashtun_se,
    non_pashtun_displaced_se,
    non_pashtun_non_displaced_se
  )),
  dv     = "econ2"
)

##### plot effects_df, Economic Status #####
effects_df$dv <- factor(
  effects_df$dv,
  levels = c("econ2"),
  labels = c("Economic Status")
)

effects_df <- effects_df[, c("dv", "group", "effect", "se")]
effects_df <- effects_df[order(effects_df$dv), ]

# Reordering the covariate and outcome categories
effects_df$group <- factor(
  effects_df$group,
  levels = c(
    "Non-Pashtuns\nNon-Displaced",
    "Non-Pashtuns\nin Displaced\nRegions",
    "Displaced\nPashtuns"
  )
)

pdf("replication/results/figures/Figure4.pdf", width = 5, height = 5)
ggplot(effects_df, aes(x = group, y = effect)) +
  geom_segment(
    aes(
      x      = group,
      xend   = group,
      y      = effect - 1.96 * se,
      yend   = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    size = 2
  ) +
  # geom_segment(aes(x = group, xend = group,
  #                  y = effect - 1.645 * se,
  #                  yend = effect + 1.645 * se,
  #                  linetype = "90% Confidence Interval")) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_wrap(~ dv, nrow = 1) +
  coord_flip() +
  xlab("") +
  ylab("") +  # "Compared to the baseline group, Non-Pashtuns Non-Displaced"
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_text(size = 16),
    strip.text         = element_text(size = 16),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom"
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
dev.off()

################################################################################
############################## Violence Exposure ###############################
################################################################################
rm(list = ls())

# Read in the exported Stata results
results <- read.csv("replication/data/intermediate/FOGHORNresults_violence.csv")
results <- sapply(results, gsub, pattern = "^=", replacement = "")
colnames(results) <- results[1, ]
results <- results[-1, ]
results <- as.data.frame(results)
results[, 2:3] <- lapply(results[, 2:3], as.numeric)

# Extract coefficients and standard errors for each group
non_pashtun_displaced_effect     <- results[results$V1 == "mean_nonpashtun_displaced_res", 2:3]
non_pashtun_displaced_se         <- results[which(results$V1 == "mean_nonpashtun_displaced_res") + 1, 2:3]
non_pashtun_non_displaced_effect <- results[results$V1 == "mean_nonpashtun_nondisplaced_res", 2:3]
non_pashtun_non_displaced_se     <- results[which(results$V1 == "mean_nonpashtun_nondisplaced_res") + 1, 2:3]
displaced_pashtun_effect         <- results[results$V1 == "mean_displaced_pashtun_res", 2:3]
displaced_pashtun_se             <- results[which(results$V1 == "mean_displaced_pashtun_res") + 1, 2:3]

# Create a dataframe with the desired groups and their effects
effects_df <- data.frame(
  group  = c(
    rep("% Respondents\nDisplaced\nPashtuns",            2),
    rep("% Respondents\nNon-Pashtuns\nin Displaced\nRegions", 2),
    rep("% Respondents\nNon-Pashtuns\nNon-Displaced",    2)
  ),
  effect = unlist(c(
    displaced_pashtun_effect,
    non_pashtun_displaced_effect,
    non_pashtun_non_displaced_effect
  )),
  se     = unlist(c(
    displaced_pashtun_se,
    non_pashtun_displaced_se,
    non_pashtun_non_displaced_se
  )),
  dv     = names(c(
    displaced_pashtun_effect,
    non_pashtun_displaced_effect,
    non_pashtun_non_displaced_effect
  ))
)

##### plot effects_df, Violence Exposure #####
effects_df$dv <- factor(
  effects_df$dv,
  levels = c("intercommunal", "sigacts"),
  labels = c("Intercommunal Violence", "Insurgent Attacks")
)

effects_df <- effects_df[, c("dv", "group", "effect", "se")]
effects_df <- effects_df[order(effects_df$dv), ]

# Reordering the covariate and outcome categories
effects_df$group <- factor(
  effects_df$group,
  levels = c(
    "% Respondents\nNon-Pashtuns\nNon-Displaced",
    "% Respondents\nNon-Pashtuns\nin Displaced\nRegions",
    "% Respondents\nDisplaced\nPashtuns"
  )
)

pdf("replication/results/figures/FigureA19.pdf", width = 7, height = 5)
ggplot(effects_df, aes(x = group, y = effect)) +
  geom_segment(
    aes(
      x      = group,
      xend   = group,
      y      = effect - 1.96 * se,
      yend   = effect + 1.96 * se,
      linetype = "95% Confidence Interval"
    ),
    size = 2
  ) +
  # geom_segment(aes(x = group, xend = group,
  #                  y = effect - 1.645 * se,
  #                  yend = effect + 1.645 * se,
  #                  linetype = "90% Confidence Interval")) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_wrap(~ dv, nrow = 1) +
  coord_flip() +
  xlab("") +
  ylab("") +  # "Compared to the baseline group, Non-Pashtuns Non-Displaced"
  theme(
    panel.background   = element_blank(),
    panel.grid.major   = element_line(colour = "gray90", linetype = "solid"),
    axis.text.x        = element_text(size = 16),
    axis.text.y        = element_text(size = 16),
    strip.text         = element_text(size = 16),
    panel.border       = element_rect(colour = "black", fill = FALSE),
    strip.background   = element_rect(colour = "black"),
    axis.title         = element_text(size = 16),
    legend.title       = element_blank(),
    legend.text        = element_text(size = 16, margin = margin(t = 20), vjust = 5),
    legend.position    = "bottom"
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
dev.off()
